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Abstract. We discuss spatial dynamics and collapse scenarios of localized waves 
governed by the nonlinear Schrodinger equation with nonlocal nonlinearity. Firstly, 
we prove that for arbitrary nonsingular attractive nonlocal nonlinear interaction in 
arbitrary dimension collapse does not occur. Then we study in detail the effect of 
singular nonlocal kernels in arbitrary dimension using both, Lyapunoff 's method and 
virial identities. We find that in the one-dimensional case, i.e. for n — 1, collapse 
cannot happen for nonlocal nonlinearity. On the other hand, for spatial dimension 
n > 2 and singular kernel ~ 1/V a , no collapse takes place if a < 2, whereas collapse 
is possible if a > 2. Self-similar solutions allow us to find an expression for the 
critical distance (or time) at which collapse should occur in the particular case of 
~ 1/r 2 kernels. Moreover, different evolution scenarios for the three dimensional 
physically relevant case of Bose-Einstein condensates are studied numerically for both, 
the groundstate soliton and higher order toroidal states with and without an additional 
local repulsive nonlinear interaction. In particular, we show that presence of local 
repulsive nonlinearity can prevent collapse in those cases. 

PACS numbers: 42.65. Jx,03.75.Kk,52.35.Mw 
1. Introduction 

Spatial and/or temporal dynamics of waves in many nonlinear systems are often 
governed by the famous nonlinear Schrodinger (NLS) equation pQ. This universal 
equation appears in many diverse physical settings including, for instance, nonlinear 
optics [2J, Bose-Einstein condensate [HI HI 15] , water waves [6] and is even discussed in the 
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context of rogue waves [7|. A particular solution to this equation - a soliton, represents 
a nonlinear localized bound state that does not change its shape when propagating or 
evolving in time and/or space [8]. The existence of solitons is brought about by the 
balance between the natural effect of dispersion or diffraction which tends to spread the 
wave and the self-focusing, provided by the nonlinear response of the medium induced 
by the wave itself. In case of optics the nonlinearity represents the self-induced refractive 
index change, while, e.g. in Bose-Einstein-Condensates (BEC) it is just the nonlinear 
potential of bosonic interaction. In both cases the presence of the nonlinearity may 
lead to self-focusing and consequently, self-trapping of light or condensate. When this 
nonlinear potential induced by the wave itself is strong enough for the self-focusing to 
dominate over the effects of dispersion/ diffract ion, the amplitude of the wave may grow 
to infinity (blow-up) within finite time or distance. This phenomenon is often referred to 
as a collapse. Collapse is a well-known effect associated with dynamics of finite waves in 
nonlinear systems. It has been studied in several physical settings including, e.g. plasma 
waves and nonlinear optics [9J, BEC [10], as well as collapse of very heavy stars into 
black holes [IT] . In the majority of works devoted to collapse the nonlinear potential 
induced by the wave is of a local character, i.e. the system response in the particular 
spatial location is determined solely by the wave amplitude (or intensity) in the same 
location. 

Here we will discuss the dynamics of waves and collapse in systems exhibiting a 
spatially nonlocal nonlinear response. In nonlocal media the nonlinear response at the 
given point depends on the wave intensity in an extended neighborhood of this point. 
The extent of this neighborhood in comparison to the spatial scale of the localized wave 
determines the degree of nonlocality [12]. It appears that spatial nonlocality is a generic 
property of many nonlinear systems and is often a result of an underlying transport 
process, such as ballistic transport of atoms [13] or atomic diffusion [TJ] in atomic vapors, 
heat transfer in plasma [15] and thermal media [T6| [T7] , or drift of electric charges as in 
photorefractive crystals [18] . It can also be induced by a long-range particle interaction 
as in nematic liquid crystals which exhibit orientational nonlocal nonlinearity [19] or 
in colloidal suspensions [20]. Recent enhanced interest in nonlocal nonlinearities has 
been stimulated by the research on Bose Einstein condensates where the inter-particle 
nonlinear interaction potential is inherently nonlocal j5j ETJ [22] . 

In contrast to local attractive nonlinearity, which leads to collapse in more than 
one dimension, nonlocal nonlinearities may actually prevent collapse. Turitsyn was the 
first to prove analytically the absence of collapse in the nonlocal NLS equation for three 
different singular nonlocal kernels in [23], including the singular gravity-like 1/r-kernel. 
In a later work Perez-Garcia et al. investigated the collapse dynamics of the nonlocal 
Schrddinger equation using variational techniques and numerical simulations [23] . They 
showed in particular, that for the nonsingular nonlocal nonlinearity collapse is always 
arrested. Absence of collapse for a singular 1/r- response and form stability of higher 
order states, such as dipole and torus, have been recently verified numerically in [25J. 
More general treatment of the NLS and ensuing dynamics and collapse is presented in 
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[25J [27]. As far as collapse in the general case of the nonlocal Schrddinger equation 
is concerned, Bang et al. showed rigorously that nonlocality arrests the collapse as 
long as the Fourier spectrum of the nonlocal nonlinear potential is absolute integrable 
[28] . Rigorous recent studies of collapse in three dimensions in case of a particular 
singular form of the nonlocal nonlinear potential have been given in [291 [301 [3Tj . Here, 
we will generalize those previous findings to arbitrary number of transverse dimensions, 
and show that for the nonlocal nonlinearity represented by singular kernels ~ l/r a , 
a = 2 there is a sharp threshold for collapse. By means of (3+l)-dimensional numerical 
simulations we elucidate different evolution scenarios, and show that an additional local 
repulsive nonlinear interaction may prevent collapse. 



2. Model 



In this paper we will consider the self-trapped evolution of the wave function tp(x,t) 
described by the nonlocal nonlinear Schrddinger equation (NNLS), given by 

idti/)(x,t) = -AVM)- (J A(x-y)\ij(y,t)\ 2 dy^(x,t), (1) 

where x, y G R n with n being the spatial dimension. A(x) is the so called nonlocal 
response or kernel which represents the nonlocal character of the nonlinearity of the 
medium. Its width determines the degree of nonlocality [12]. In particular, for A(x) ~ 
5(x) the above equation describes the standard Kerr-type local medium. While this is 
only a phenomenological model it nevertheless describes very well general properties of 
nonlocal media. Moreover, for certain forms of the nonlocal response function this model 
represents actual experimentally realizable physical systems. This is the case of nematic 
liquid crystals [HJ [32] or thermal nonlinearity in liquids [17J and plasma [15] where, in 
one dimensional case A(x) ~ exp(— \x\) and in two dimensions A(x) is proportional to 
the modified Bessel function of the second kind [33]. In addition, Eq. ([[]) (also known 
as Gross-Pitaevskii equation) governs the spatio-temporal dynamics of Bose-Einstein 
condensates where the analytical form of the nonlocal Kernel is a consequence of the 
dipolar character of the long-range particle-particle interaction [221 EH ESI [36J. While, 
in general, the nonlocal interaction potential in condensates is of rather complex and 
anisotropic character, under certain experimental conditions it may acquire "gravity- 
like" 1/r form [3?]. 



3. Analytical estimates for possible evolution scenarios 

In a vivid picture, collapse refers to a phenomenon when the amplitude of the wave 
function ip goes to infinity within finite time or distance (for stationary processes). 
However, in the mathematical context collapse is usually defined by means of a diverging 
gradient norm || V^>||2 (see Appendix A for notations). For practical considerations, this 
is essentially the same, because a wave function maintaining a finite L 2 norm, \\1pW2, 
attains infinite amplitude when its gradient norm || V?/' || 2 diverges [38] . 
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3.1. Lyapunoff methods 

In this subsection, we will use Lyapunoff 's method [39J, to discuss collapse scenarios. 
The idea is, roughly speaking, to show that the Hamiltonian 

H = IS \ V ^ dxdy -\ II \i>(x)\ 2 A(x-y)\ij(y)\ 2 dxdy (2) 

for Eq. (0Q) is bounded from below. Then, since in our case the L 2 norm \\1pW2 arid the 
Hamiltonian % are conserved quantities, it is possible to show that the gradient norm 
||V^||2 is bounded from above, which implies that it cannot diverge. Moreover, we 
can expect the existence of at least one stable soliton solution when the Hamiltonian is 
bounded (see, e.g., [40J). 

We can decompose a general kernel A = A + + A", A + > and A~ < 0. Then, the 
kernels A + and A~ represent attraction (in BEC) or focusing (in nonlinear optics) and 
repulsion or defocussing, respectively. Since the amplitude of the wave function has to 
become infinity at a certain point, one can expect that a singular kernel A + promotes 
collapse. Using Young's convolution theorem, one can immediately infer the following 
related fact. Without loss of generality, we restrict our following considerations to solely 
attractive kernels A = A + . 

Lemma 1. Let A G L°°(R n ), ij){t = 0,x) G /f 1 (R n ). Then no collapse occurs at all 
times, since HV^Hi ^ s bounded from above. 

Proof. 

«>l|vvil3-^iHlMoo. (3) 

Here, C = 1 is a sharp constant (see, e.g., [H|). Since HV'lh and the Hamilton % 
itself are conserved quantities under the flow of time (or propagation), one can read the 
inequality ([3]) as bounding HV^H! from above in terms of conserved quantities. □ 

Let us try to connect the inequality with the earlier results of the stability of 
finite beams governed by the NNLS equation (JTJ) [28J. This work states that in order to 
arrest collapse the spectrum of the nonlocal kernel A should be absolutely integrable. 



Since ||^4||oo < C||^l|i> where A denotes the Fourier transform of A (see Appendix A ), 
the inequality ([3]) clearly contains the result of Ref. [28] and also explains earlier reports 
on collapse arrest for nonsingular kernels observed in numerical simulations [21]. 

Further on, inequality fl3]) shows that collapse can only be relevant if the kernel 
A diverges at a certain point (usually at the origin x = 0). Then the assumption 
A G L°°(R n ) is violated, as is the case, e.g., of local interaction [A ~ S(x)]. However, a 
diverging kernel A does not necessarily imply collapse. In order to ease the restrictions 
on the nonlocal kernel which ensure the collapse arrest we have to restore to the gradient 
norm ||V^||2 apart from the L 2 -norm ||^||2- As we will show below one can weaken the 
sufficient restriction A G L°° by employing the Hardy-Littlewood-Sobolev inequality 
(see |Appendix B[ ) and the Gagliardo-Nirenberg-Sobolev inequality (see |Appendix Cp . 
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Lemma 2. Let x G R n , A = f(x)/r a , r = \x\ with < f(x) < oo, a < min(2,n), 
ip(t — 0, x) G H 1 (R, n ) ! K := max x f < oo. Then collapse cannot occur. 



Proof. 



^ > IIV^HI - -K 



M*)IXy)l ; 

\x — y\ a 



■dydx 



> mm - x -Kcut 2n/{n _ a/2) = - \*c ( / m 2n/{n - a/2) dx 



,/2 



(4) 



In the second step we just applied the Hardy-Littlewood-Sobolev inequality ( IB.lj) . 
Introducing o = a/(2n — a), we can rewrite the last expression in (j3J) and find: 



n>\\v^\\l--KC 



>\\V1>\\I--kC 



2 

-r 



K 



2n/(n-o/2) dx 

|V^| 2 da; 



,/2 



-kC 



\2+2a 



dx 



2 



dx 



2+tr(2-n) 
2 



\m\\i-^c\\vmmtr 



(5) 



Here, we employed the Gagliardo-Nirenberg-Sobolev inequality (IC.lj) . and C := 
KCK OL l na is some positive constant independent of ip. Since a < 2 and because the 
Hamiltonian H and the L 2 -norm HV'lh are conserved quantities, inequality (jSJ) can be 
read as a uniform bound of the gradient norm || 2- D 

For dimension n — 1, the convolution term in the Hamiltonian 7{ and expression (j3J) 
diverge without regularization for singular kernels with a > 1, and hence singularity 
cannot be strong enough to promote collapse. Generally, when regularizing the kernel 
by, e.g., introducing a cut-off, again the wave function does not collapse since the kernel 
is then nonsingular. 

If we apply the lower bound for the Hamiltonian (jSJ) for the particular exponent 
a = 2 we find 

n>\m\i(i-lcu\\i). (6) 



Because we know the sharp constant C := KCK a / n(T [k = 1, C = C(n,a) according 
to Eq. (IB.2p . K = K(a,n) according to Eq. (10. 2ft ] we can use the inequality (JHD 
to compute a lower bound for the "critical mass" \\1pW2 = 2/C for collapse. In the 
physical relevant case of n = 3, we find this lower bound numerically as 2/C = 2.80 0, 
which is slightly lower than the numerical result for the mass of the groundstate soliton, 
HV'Gslli = 2-85 (see Sec. 14. ip . This finding may seem surprising, because for the local 



I Note that the constant K 
numerically. 



K(a,n) involves the quantity ||-R||2 CT , which has to be computed 
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nonlinear Schrddinger equation the exact soliton mass can be computed from sharp 
constants in the lower bound of the Hamiltonian H21 03]. However, a closer look at 



the Hardy-Littlewood-Sobolev inequality (Sec. Appendix B ) reveals that equality can 
be expected only for wavefunctions ip of the form (IB. 31) . Since our groundstate solitons 
t/>gsOe) are not of the form (IB. 31) I§1. we have to expect 2/(7 < H^gsIII- Finally, we verified 
numerically that for the groundstate soliton (IB.lj) is indeed an inequality. Interestingly, 
then the ratio of r.h.s. and l.h.s. of (lB~Tj) is about 2.85/2.80 = ||V>gs||! x C/2. 



3.2. Virial identities 

The idea to use virial identities to derive sufficient conditions for collapse was first 
established in the context of the local two-dimensional nonlinear Schrddinger equation 
by Vlasov, Petrishchev, and Talanov [H], and the resulting condition H < and is 
often referred to as the VPT criterion. The original idea behind the VPT criterion has 
been extended to various modifications of the original NLS equation, e.g., the Gross- 
Pitaevskii equation for a BEC with dipole-dipole interactions [29] or, more general, the 
NLS with l/r a nonlocal response in three dimensions [31]. Following the calculations 
presented in |31j . one may observe that the derived virial identities hold for arbitrary 
dimension. The virial is given by V(t) := J r 2 \%l)\ 2 dx with r = \x\. 

Lemma 3. Consider a kernel A = l/r a , with x 6 R n and a < n. Then the second time 
derivative of the virial in arbitrary dimension n is given by 



^- 8 (/lv^-j//^) 



(7) 



\x - y\ 



Proof. This task was already done for three dimensions in [31]. Therefore we only briefly 
sketch the derivation of this virial identity here. Apparently, 

d t V = 43 J ij*xVipdx (8) 

Using the equation of motion, a longer but straightforward calculation leads to (see also 

d tt V = 8 [ |V^| 2 -4n/ l ^ )|2| ^ )|2 dydx + 4 /*V(hK*)| 2 ) f P^dydx (9) 
J J \x-y\ a J J \x-y\ a 

§ lim ipGs ~ exp(— r) r>c r~ 4 , with |x| = r. 
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Integrating the last term by parts and re-writing it in a more favorable way gives 

iV(|</>W| 2 ) j ^^-Ayix (10) 

= „ ft KWMfJL - H // ,(iqM£)JW dsd:c (n) 
J J \x y\ ^ J J \*c y\ 

2 yy y |x-i/i-+ 2 y 

Plugging this into Eq. fl9]), the first term cancels out, and the two symmetric terms with 
respect to x and y add up and give the desired expression. □ 

In particular, Lemma |3] implies duV = 8H for A = 1/r 2 (or a = 2) in arbitrary 
dimension. Together with a particular form of Holder's inequality 

r\ip\\Vip\dx < r*\i()\ 2 dxj (^J |V^| 2 dx^ , (12) 
we can derive a sufficient criterion for collapse. 

Lemma 4. Let x G R n , n > 2, A = 1/r 2 , ip(t = 0,x) G if 1 (R n ), and the Hamiltonian 
7i < 0. T/ien i/ie solution ip undergoes collapse. 

Proof. As seen above, the second time derivative of the virial V equals eight times the 
Hamiltonian. The virial is positive for t = and has the apparent meaning of the mean 
width of ip. Hence, for a negative Hamiltonian, the virial vanishes at a certain time t. 

Using partial integration and Holder's inequality f|T2l . one can infer that a vanishing 
virial leads to collapse due to conservation of the L 2 -norm \\1pW2 of the wavef unction and 

J \^\ 2 dx<C^J r 2 \^\ 2 dx^j ' (J iWfdx^) ' , (13) 

where the positive constant C does not depend on ip. □ 

In other words, if the mean width or the virial becomes zero, the gradient norm 
has to diverge in order to keep the L 2 -norm constant. Closer inspection of Lemma [3] 
reveals that a negative Hamiltonian also implies collapse for a > 2, as already pointed 
out in [31] . Together with Lemma |2l we see that a = 2 is a sharp boundary for collapse 
in arbitrary dimension n > 2. 

3.3. Self-similar solutions for the 1/r 2 kernel 

Let us now have a look at the special case A = 1/r 2 . As we have seen in the previous 
section, this singular kernel marks the boundary between collapsing and non-collapsing 
solutions. This is somewhat similar to the two-dimensional local nonlinear Schrddinger 
equation, and we will find many aspects of this famous equation here. However, in the 
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nonlocal case the collapse threshold does not appear as a critical dimension, but as a 
critical exponent of the singularity of the kernel. 

The NNLS equation with kernel A = 1/r 2 appears as equation of motion of the 
following Lagrangian density, 

C:= l -^d t r-rd^) + \^^\ 2 - l -W\ 2 j j^^\^(y,t)\ 2 dy. (14) 

To simplify the considerations and following j9], let us assume radial symmetry. A 
possible ansatz for a self-similar solution (or similariton jl5]) is 

V> = B(t)R(r/a(t))e i6 ^ (15) 

with r = |x|. Here, R(t/a(t)) is a real profile, and a(t) is the wave radius. Due to the 
conservation of the L 2 -norm \\ip\\2 and the related continuum equation, which comes 
about from the global U(l) invariance of the Lagrangian, one finds the interdependence 
between B and a: 

\B\ 2 a n =: P = const. (16) 
An appropriate ansatz for the radially symmetric phase is 

9(r,t)=9 2 (t)r 2 + 9 (t). (17) 

Using 9 and P, one may express the spatially integrated Lagrangian density as follows: 



L = -\B\ 2 a n 
where we introduced 



Kod t e Q + K 2 a 2 (dt0 2 + A9 2 2 ) + -L - (5 \ B \ 4 a 



\i„2n-2 



K m := J r m R 2 (r)dx, (3 := 1/2 J ^M^^- dxdy, and 7 := J (d r R(r)) 2 dx. (19) 

The variational problem 5L = with respect to B, , # 2 and a gives the equations of 
motion for these parameters. Using P = \B\ 2 a n and 9 2 = (dta)/4a, one may rewrite the 
equation of motion for a in the form 

d u a = r r. (20) 

Elementary integration with the initial conditions Oq := a(t = 0), and d t a(t = 0) = 
gives 

a(t) = a Jt 2 (^ 1 -^)+l. (21) 
Now one can estimate the critical time t cr , where a — > as 
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In terms of P, one may express the condition for collapse as P > %. Compared to 
the virial results from the previous section, this threshold value indicates that the 
Hamiltonian H becomes negative. 

It is important to stress here that the variational approach presented above strongly 
relies on the ansatz Eqs. ( fl5l) and (fTTl) . A priori, there is no guarantee that variational 
results show even qualitative agreement with the exact solutions. A classical example 
of unsuccessful application of the variational method is given in [46J, where the authors 
show that the description of solitons interacting with radiation in the one dimensional 
local NLS equation fails. However, as far as collapse in the two and three dimensional 
local NLS is concerned, the above ansatz is capable of giving useful estimates [4T1 l4~8l l9]. 
As we will see in the following sections, in our case variational predictions show some 
qualitative agreement with rigorous numerical simulations in the physically relevant case 
of n = 3. In the light of Ref. [IE], this is indeed the best we can expect. 

4. (3+l)-dimensional numerical investigation of collapse dynamics and 
collapse arrest for the 1/r 2 kernel 

In this section we will investigate the dynamics of finite wave packets using numerical 
integration of the NNLS equation ([1]). We start by using variational methods to find 
families of approximate localized solutions to the stationary version of this equation, 
and compare them to their numerically exact counterparts computed by imaginary time 
evolution [19] or iteration technique [33J. Both of them will be subsequently used as 
initial conditions in numerical simulations in order to illustrate possible scenarios of 
time evolution of the wave function in case of the nonlocal 1/r 2 kernel. Our (3+1)- 
dimensional simulation code is based on a Fourier split-step method. In order to resolve 
the singularity of the kernel properly, which is of crucial importance when studying 



collapse phenomena, we use a method outlined in Appendix D In the following, we 
will also investigate the impact of an additional local repulsive (defocusing) term on 
the collapse dynamics. In particular, we will be interested whether the presence of the 
repulsive local interaction can arrest the collapse. To this end, let us consider a NNLS 
equation of the following form: 



id t ip(x,t) = —Aip(x,t) — 



' ^(y,t)\ 2 dy-5mx,t)\ 2 



|2 



^(x,t). (23) 



\x - y\ 

where the parameter 5 in front of the local term is either or 1. On physical grounds, 
the additional local repulsive term occurs in Bose-Einstein-Condensates as a result of 
s-wave scattering [5J. 

4-1. Approximate variational and exact numerical solutions for the groundstate 

To find approximate localized groundstate solutions to the NNLS equation we use the 
following ansatz 



1)(x,t) := A exp [ -— ) exp(tEt) (24) 
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Figure 1. Existence curves of the groundstate soliton [(a) and (b)] obtained by 
variational analysis and numerical computation for nonlocal attractive 1/r 2 response 
(blue solid and dashed green line, resp.) and with competing local repulsive interaction 
(black solid and dotted blue line, resp.). c) Time evolution of the amplitude of the 
numerically exact groundstate. If the mass of the groundstate is increased, the solution 
collapses (shown here for yLOl and 1.01 x V'gs), an d if the mass is decreased, the 
solution spreads upon propagation (shown here for 0.99 and a/0.99 x i/'gs)- The inset 
shows a zoom into the early propagation stage. 

and vary the corresponding Lagrangian with respect to width o and amplitude A . Here, 
we restrict our considerations to dimension n = 3, so that x = (xi,X2,Xs). The results 
are shown in Figs. H^a-b), which depict the so called existence curves for the solutions, 
i.e., the relations between mass M := \\1pW2, energy E and width a. 

The solid lines in Fig. DJa) represent families of solutions obtained by variational 
analysis with (black line) and without (blue line) repulsive local interaction (5 = 0), 
and the dashed green horizontal line together with the blue dots show the same for 
the numerical exact soliton solution of Eq. ( 123]) . We expect from this plot, that the 
groundstate may indeed collapse without repulsion (5 = 0), since for fixed mass M the 
width a can be arbitrarily small. Typical length-, time- and amplitude-scales R c , T c 
and ip c for an exact solution behave like 1/T C ~ 1/i?^ ~ ^ 2 c Rc- Hence, the invariance of 
the mass M with respect to energy E can be expected from M ~ ipc-^-c ~ const. The 
same argument holds for arbitrary dimension n, where we have 1/R 2 ~ ip^R^/R^. In 
contrast, when including the additional local repulsive interaction competing with the 
nonlocal attractive one, we expect absence of collapse from Fig. [Ha), since the mass M 
of the groundstate is increasing with decreasing width a. 

Before presenting rigorous numerical simulations of the time evolution of the above 
solutions, let us evaluate the expression for the critical collapse time t cr derived in 
section 13.31 Plugging the Gaussian approximate variational groundstate as well as the 
numerically exact soliton as self-similar profile into Eq. (122]) . the critical time t cr diverges, 
and we expect that collapse happens at infinite times. However, when multiplying the 
profiles by some factor larger than one, t cr becomes finite, and we expect collapse at 
time t cr . 

These predictions are confirmed by our numerics for the exact groundstate, as can 
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Figure 2. The evolution of the amplitude for the initial variational groundstate 
is shown in a). The blue circles correspond to snapshots of the evolution shown in 
Fig. [3] In b), the evolution of the variational groundstate multiplied by 2 (left curve) 
and 1.5 (right curve) is plotted and compared to the estimates for the critical time 
t CI . In c), the possibility to prevent collapse using a local repulsive term is illustrated: 
The variational groundstate multiplied by a factor of 5 remains much smaller than the 
collapsing ones without repulsion shown in b). 

be seen in Fig. H^c). The groundstate is self-trapped and evolves with small change 
in amplitude (invisible in the figure) due to numerical imperfections. If we multiply 
the exact groundstate by a factor smaller than one (e.g. V0.99), it spreads with time, 
whereas if we multiply it by a factor larger than one (e.g. VL01), it collapses. The 
critical times t cr obtained from Eq. (122]) are in excellent agreement with the simulations. 
In contrast, using the approximate variational groundstate as initial condition, we see 
that it collapses [Fig. (2K)] after rather short time. Thus, the estimate for the critical 
time t CT = oo is incorrect. This is related to the fact that the approximate variational 
Gaussian profile has a larger mass than the exact soliton (M = 3 > Mqs)- When 
we multiply the variational groundstate by a factor larger than one [see Fig. 121(b)], we 
see collapse earlier than predicted by Eq. ( 122]) . Collapse-times t CI are systematically 
overestimated. Finally, collapse can be prevented by adding a local repulsive term to 
the nonlinear potential in the original NNLS equation [see Fig. 121(c)]. 

The time evolution of the initially Gaussian variational groundstate in Fig. Ufa) 
can provide us with some deeper insight into the collapse dynamics. From our previous 
attempts to predict the critical collapse time t CT , we may suspect that the exact soliton 
itself is the self-similar profile. In order to corroborate this hypothesis, Figs. [3fa), (b) 
and (c) show some snapshots of the radial profiles during time evolution corresponding 
to Fig. 121(a). Already after evolving over a short time, the wave function clearly deviates 
from a Gaussian profile, and resembles more the numerically exact solution [see black 
curves in Fig. 121(b)]. Further propagation shows that the collapsing solution converges to 
the exact groundstate. A similar behavior is know for the classical two-dimensional local 
nonlinear Schrddinger equation, where the famous Townes soliton naturally emerges 
in the collapse process (see, e.g., [50l El]). Moreover, we checked that an initially 
ellipsoidal wave function with slightly too large mass will collapse as a radially symmetric 
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x i x i x i 

Figure 3. Snapshots of the time evolution of the initially Gaussian variational 
groundstate of Fig. [^a) at a) t = 0, b) t = 1.95 and c) t = 2.25. The black curves show 
numerically exact solitons with the same amplitude. Clearly, some part of the initially 
Gaussian wave function is radiated away, and the actually collapsing part converges 
towards the exact soliton. 



formation, that again resembles the numerically exact solution (not shown). 
4-2. Evolution scenarios for toroidal states 

Finally, we want to have a look at higher order solitons and their collapse behavior. 
For this purpose, we consider toroidal states, i.e. the three-dimensional generalization 
of a vortex. Toroidal states have an azimuthally symmetric density together with an 
azimuthal phase ramp, which leads to a phase singularity of integer topological charge 
m in the origin 

ip(r, <f>, z, t) = p(r, z) exp(im0) exp(iEt). (25) 

In order to find numerically exact toroidal solitons (or vortices), we use imaginary time 
evolution starting from tori computed from a variational ansatz. As for the groundstate 
soliton in Sec. I4.1[ we find that adding 1% to the mass M leads to collapse of the 
toroidal soliton, whereas subtracting 1% leads to derealization of the wavefunction, 
that cannot self-trap anymore [see Fig. IU(a)]. However, the toroidal state is azimuthally 
unstable. Adding a small amplitude noise is enough to trigger the instability and due 
to the interplay between azimuthal phase ramp and amplitude modulation along the 
ring the whole structure starts to rotate [52]. Upon further propagation, the ring 
breaks into two humps (groundstate solitons), that eventually collapse individually [see 
Figs. BJb)-(f)]. This secondary collapse is possible because the mass of the torus fulfills 
M Torus > 2 x M GS . 

5. Conclusions 

The phenomenon of collapse in the nonlocal nonlinear Schrddinger equation has been 
investigated in arbitrary dimensions. We showed that nonsingular nonlocal kernels do 
not support collapsing solutions. More general, it has been shown that for singular 
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Figure 4. Time evolution scenarios for the single charged (m — 1) toroidal soliton. 
a) Time evolution of the amplitude. When reducing the mass by 1%, the solution 
spreads, whereas it collapses for masses M > 7.98. b)-f) The vortex may decay into two 
groundstate solitons due to an azimuthal instability, which then collapse individually, 
g) The azimuthal instability can be triggered by adding random noise (here 1% in 
amplitude), so that instead of spreading (blue solid line) a torus with M = 7.90 decays 
and collapses (black dashed line). The corresponding 3D evolution of the modulus 
square of the wave function in the (xi,X2) plane (2:3 = 0) is depicted in b)-f). 

kernels ~ l/r a , the necessary condition for collapse to occur is a > 2. Various collapse 
scenarios for the 1/r 2 kernel have been studied numerically in the physically relevant case 
of three transverse dimensions (n = 3). The critical time (or distance) when collapse 
occurs has been estimated by using self-similar solutions. Apart from the groundstate 
solution, collapse of a toroidal state with absence and presence of local repulsive 
interaction has also been studied in this particular system. Generally, it appears that 
the initial torus is azimuthally unstable. Finally, we showed that independently of the 
initial state, collapse can be arrested by adding a local repulsive contribution to the 
attractive 1/r 2 kernel. 
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Appendix A. Notations 

We use common functional analysis definitions. The symbol L p (R n ) denotes the space 
of measurable functions for which the Lebesgue measure 




(A.l) 
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exists. A function is in L^R"), if there is a constant C > 0, such that \if)\ is smaller 
than C for almost every x G R n . The symbol H 1 (R n ) (Sobolev space) is the space of 
measurable functions for which the norm, defined as 



H 1 



IIWIh + IHh 



(A.2) 



is smaller than oo. Here, J denotes 

The Fourier transform / of a function / is given by 



"OO 

-co' 



f(k) = J e- 2i7T{k ' x) f{x)dx 

with x,k G R n , and (•, •) denotes the scalar product. 

Apart from that, \x\ = \fx\ + x\ + . . . = r is the usual Euclidean distance. 



(A.3) 



Appendix B. Hardy-Littlewood-Sobolev inequality 

Let p, I > 1 and < a < n with 1/p + a/n +1/1 = 2. Let f e L p (R n ) and h G L ; (R"). 
Then there exists a sharp constant C(n,a,p) independent of / and h, such that 

f(x)h(y) 




\x - y\ 



-dxdy 



<C(n } a,p)\\f\\ p \\h\\i- 



(B.l) 



The explicit functional dependence of the C(n,a,p) is known (see e.g. [H]). For our 
purposes, the special case p = I = 2n/(2n — a) is important, where the sharp constant 
is given by 

~T(n/2y- 1+a/n 



C{n,a,p) =C(n,a) = n a/2 



T(n/2 -a/2) 



Tin) 



T(n - a/2) 

In this case there is equality in (IB. ip if and only if h = Bf and 

f(x) = A{ 1 2 + \x-a\ 2 y (2n - a)/2 . 
for some complex constants A and B, ^ 7 real and a G R™. [H 



(B.2) 



(B.3) 



Appendix C. Gagliardo-Nirenberg-Sobolev inequality 

Let / G H 1 (R, n ) and < a < 2/(n — 2). Then there exists a sharp constant K(a,n) 
independent of /, such that 

an 2+g-(2-ra) 

\f(x)\ 2 ° +2 dx < K(a,n) ( f |V/(a;)| 2 dx^ ' ( [ \f(x)\ 2 dx] ' . (C.l) 



The explicit functional dependence of the K(a,n) is known (see e.g. jlj), and reads 



K(a,n) = {<t + L 



2 (2 + 2a -an 



-l+cro/2 



[an 



an/2 



\ K \\2 



(C.2) 



where R is the positive solution of AR — R + R 2a+1 = 0. 
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Appendix D. Numerical implementation of the singularity of 1/r 2 

Since we study collapse phenomena, it is crucial to define the properly in a numerical 
sense kernel 1/r 2 and, in particular, its value at r = . To this end we used the following 
procedure (see, e.g., [53]). Firstly, we decompose the singular kernel into a short range, 
singular at the origin, and a nonsingular, long-range contributions. 



The long-range part can be used directly as it is in a real space, whereas the short-range 
part ( e ri ° ) will be treated in the Fourier domain. In our numerics, the parameter wq 
was chosen to be Wq = 7Ax/2tc, with Ax the step size of the spatial grid. In Fourier 
domain, the product of two functions becomes a convolution, that is in particular also 
well defined for certain singular functions such as l/r a with a < n, if the other function 
is "well behaved". Hence, in Fourier domain the short ranged contribution can be easily 
calculated, and the limit \k\ — > is then well defined. 
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